##################################
# Map of Chile, all visas 
##################################

rm(list=ls())

# packages
require(maptools)
require(descr)
require(RColorBrewer)
require(plotrix)
library(foreign)
library(tibble)
library(ggplot2)

#sink("20_map_all.txt")

# Load shapefiles
divcomunas <- readShapePoly(fn='final_maps_2023april6.shp')
names(divcomunas)
#View(divcomunas@data)

###############################
# prepare plot
###############################

# ordering the vector by exposure indicator
divcomunas@data <- divcomunas@data[order(divcomunas@data$t_stagg_a0),]

# colors based on exposure
divcomunas@data$colors = 0
divcomunas@data$colors[divcomunas@data$t_stagg_a0==0] = "royalblue1"
divcomunas@data$colors[divcomunas@data$t_stagg_a0==1] = "indianred1"
divcomunas@data$colors[divcomunas@data$t_stagg_a0==2] = "lightgrey"
table(divcomunas@data$colors)

# going back to the original order of the data in the shapefile
divcomunas@data <- divcomunas@data[order(divcomunas@data$order),]
head(divcomunas@data)

# save plot
tiff('figureA3.tiff', units="in", width=10, height=7, res=600, compression = 'lzw')
plot(divcomunas, col=as.character(divcomunas@data$colors), lty=0)
legend(x="bottomleft",
       legend=c("Control", "Exposed", "Not included in CEP"), 
       fill = c("royalblue1","indianred1","lightgrey"))
dev.off()

sink()

